Priors
\(\alpha\)
tibble(x = seq(0, 1.3, length = 10^(5)),
y = with(prior_params, gamma_density(x,
mean = alpha_mean,
sd = alpha_sd,
bounds = alpha_bounds))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $\\alpha$"),
fill ='',
subtitle = paste0("Mean: ", prior_params$alpha_mean,
", SD: ", prior_params$alpha_sd))\(\beta\)
tibble(x = seq(0, 1, length = 10^(5)),
y = with(prior_params, beta_density(x,
mean = beta_mean,
sd = beta_sd,
bounds = beta_bounds))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $\\beta$"),
fill ='',
subtitle = paste0("Mean: ", prior_params$beta_mean,
", SD: ", prior_params$beta_sd))\(P(S_1|\text{untested})\)
tibble(x = seq(0, 1, length = 10^(5)),
y = with(prior_params, beta_density(x,
mean = s_untested_mean,
sd = s_untested_sd,
bounds = s_untested_bounds))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $P(S_1|untested)$"),
fill ='',
subtitle = paste0("Mean: ", prior_params$s_untested_mean,
", SD: ", prior_params$s_untested_sd))\(P(S_0| \text{test}_+, \text{untested})\)
tibble(x = seq(0, 1, length = 10^(5)),
y = with(prior_params, beta_density(x,
mean = p_s0_pos_mean,
sd = p_s0_pos_sd,
bounds = p_s0_pos_bounds))) %>%
ggplot(aes(x=x, y = y)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0,ymax=y),
fill="black",
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $P(S_0|test_+,untested)$"),
fill ='',
subtitle = paste0("Mean: ", prior_params$p_s0_pos_mean,
", SD: ", prior_params$p_s0_pos_sd))COVID-19 Trends and Impact Survey
ctis_smoothed <- tar_read(ctis_smoothed,store = here("_targets"))\(\beta\)
ctis_smoothed %>%
filter(keep) %>%
mutate(state=toupper(state)) %>%
select(date,
state,
imputed_beta,
beta_estimate_smoothed,
beta_estimate_spline_smoothed) %>%
pivot_longer(contains("beta")) %>%
mutate(name = case_when(
name == "beta_estimate_smoothed" ~ "LOESS smoothed",
name == "beta_estimate_spline_smoothed" ~ "Spline smoothed",
name == "imputed_beta" ~ "Survey Value"
)) %>%
ggplot(aes(x=date, y=value, color = name,
alpha = name, linewidth=name)) +
geom_line() +
facet_wrap(~state, ncol=4) +
scale_alpha_manual(values=c("LOESS smoothed" = .9,
"Spline smoothed" = .9,
"Survey Value"=.3),
name='') +
scale_linewidth_manual(values = c("LOESS smoothed" = 1.05,
"Spline smoothed" = 1.05,
"Survey Value"=.5)) +
scale_color_manual(values=c("#3381FF", "#B58746", "#26900F"),
name ='') +
guides(alpha="none",
linewidth="none",
color = guide_legend(override.aes = list(linewidth = 3,
alpha = c("LOESS smoothed" = .9,
"Spline smoothed" = .9,
"Survey Value"=.3)),
nrow=3)) +
theme_c(legend.position="top") +
ylim(0,1) +
scale_x_date(date_breaks="3 months",
date_labels = "%b %Y") +
labs(y = TeX("Survey Estimate of $\\beta$"),
x= "",
title = TeX("Comparing Approaches for Smoothing Survey Estimates of $\\beta$"))The ratio of the screening test positivity over the overall test positivity from the COVID-19 Trends and Impact Survey is taken to be the estimate of \(\beta\). We compare two approaches to smoothing: cubic spline smoothing with 2 knots (July 15th, 2021 and December 1st, 2021) to LOESS smoothing with a span of 0.33.
\(\Pr(S_1 | \text{untested})\)
ctis_smoothed %>%
mutate(state=toupper(state)) %>%
select(date,
state,
contains("s_untested")) %>%
pivot_longer(contains("s_untested")) %>%
mutate(name = case_when(
name == "s_untested_smoothed" ~ "LOESS smoothed",
name == "imputed_s_untested" ~ "Survey Value"
)) %>%
ggplot(aes(x=date, y=value, color = name,
alpha = name, linewidth=name)) +
geom_line() +
facet_wrap(~state, ncol=4) +
scale_alpha_manual(values=c("LOESS smoothed" = .9,
"Survey Value"=.6),
name='') +
scale_linewidth_manual(values = c("LOESS smoothed" = 1.02,
"Survey Value"=.9)) +
scale_color_manual(values=c("Survey Value" ="black", "LOESS smoothed"="darkred"),
name ='') +
guides(alpha="none",
linewidth="none",
color = guide_legend(override.aes = list(
linewidth = 3,
alpha = c("LOESS smoothed" = .9,
"Survey Value"=.3)),
nrow=3)) +
theme_c(legend.position="top") +
scale_x_date(date_breaks="3 months",
date_labels = "%b %Y") +
labs(y = TeX("Survey Estimate of $Pr(S_1|untested)$"),
x= "",
title = TeX("Comparing Approaches for Smoothing Survey Estimates of $Pr(S_1|untested)$"))The percentage of the population experiencing COVID-19-like illness is taken to be the estimate of $(S_1|). The LOESS smoothed estimate with a span of 0.2 is shown in red.
Distributions of \(\beta\) by State
# raw ctis data
ctis_raw <- readRDS(here('data/data_raw/ctis_all_states.RDS'))
states_keep <- ctis_smoothed %>%
filter(keep) %>% pull(state) %>% unique()
ctis <- ctis_raw %>%
select(signal, state=geo_value, date, value) %>%
pivot_wider(names_from=signal, values_from=value) %>%
filter(state %in% states_keep) %>%
mutate(beta=smoothed_wscreening_tested_positive_14d/
smoothed_wtested_positive_14d,
state=toupper(state))ctis %>%
group_by(state) %>%
summarize(median = median(beta, na.rm=TRUE),
mean = mean(beta, na.rm=TRUE),
q1 = quantile(beta, .025, na.rm=TRUE),
q2 = quantile(beta, .975, na.rm=TRUE)) %>%
ggplot(aes(x= fct_reorder(state, mean, .desc=TRUE),
y = mean)) +
geom_errorbar(aes(ymin=q1,ymax=q2),width=.2) +
geom_point(color="darkred", size=2) +
scale_y_continuous(n.breaks=10) +
labs(y= TeX("Empirical Estimates of $\\beta$"),
x="State",
title = TeX("Empirical Estimates of $\\beta$ from the COVID-19 Trends and Impact Survey"),
subtitle = "2.7% Percentile, 97.5% Percentile, and Median By State" ) +
theme_c(axis.title=element_text(size=12))Since \(\beta\) represents the ratio of \(\Pr(\text{test}_+|\text{untested}, S_0)\) to \(\Pr(\text{test}_+|\text{tested})\), we estimate \(\beta\) from the COVID-19 Trends and Impact Survey as the ratio of the screening test positivity over the overall test positivity. Here, we consider the 2.5% and 97.5% percentiles and mean for each state. Most have a mean near 0.2.
ggsave(here("figures/emp-estimate-beta-dist-by-state.pdf"))#
# Logistic transformation of the Beta CDF.
#
f.beta <- function(alpha, beta, x, lower=0, upper=1) {
p <- pbeta(x, alpha, beta)
log(p/(1-p))
}
#
# Sums of squares.
#
delta <- function(fit, actual) sum((fit-actual)^2)
#
# The objective function handles the transformed parameters `theta` and
# uses `f.beta` and `delta` to fit the values and measure their discrepancies.
#
objective <- function(theta, x, prob, ...) {
ab <- exp(theta) # Parameters are the *logs* of alpha and beta
fit <- f.beta(ab[1], ab[2], x, ...)
return (delta(fit, prob))
}
#
# Solve two problems.
#
par(mfrow=c(1,2))
alpha <- 15; beta <- 22 # The true parameters
for (x in list(c(1e-3, 2e-3), c(1/3, 2/3))) {
x.p <- f.beta(alpha, beta, x) # The correct values of the p_i
start <- log(c(1e1, 1e1)) # A good guess is useful here
sol <- nlm(objective, start, x=x, prob=x.p, lower=0, upper=1,
typsize=c(1,1), fscale=1e-12, gradtol=1e-12)
parms <- exp(sol$estimate) # Estimates of alpha and beta
#
# Display the actual and estimated values.
#
print(rbind(Actual=c(alpha=alpha, beta=beta), Fit=parms))
#
# Plot the true and estimated CDFs.
#
curve(pbeta(x, alpha, beta), 0, 1, n=1001, lwd=2)
curve(pbeta(x, parms[1], parms[2]), n=1001, add=TRUE, col="Red")
points(x, pbeta(x, alpha, beta))
}## alpha beta
## Actual 15.00000 22.00000
## Fit 14.99983 21.99811
## alpha beta
## Actual 15.00000 22.00000
## Fit 14.99926 21.99926
ctis %>%
group_by(state) %>%
summarize(median = median(beta, na.rm=TRUE),
mean = mean(beta, na.rm=TRUE),
q1 = quantile(beta, .025, na.rm=TRUE),
q2 = quantile(beta, .975, na.rm=TRUE)) %>%
ggplot(aes(y=fct_reorder(state, mean, .desc=TRUE))) +
ggridges::geom_density_ridges(
aes(x=beta,
y=fct_reorder(state, beta, .fun =median)),
data=ctis,
fill="#596281",
quantile_lines=TRUE,
quantiles = c(.025,.975),
quantile_fun=function(beta,...) c(
# quantile(beta,.025, na.rm=TRUE),
median(beta,na.rm=TRUE)
# quantile(beta,.975, na.rm=TRUE)
)) +
# geom_errorbar(aes( xmin=q1,xmax=q2),width=.2) +
scale_x_continuous(n.breaks=10, limits=c(0,1)) +
labs(x= TeX("Empirical Estimates of $\\beta$"),
y="State",
title = TeX("Empirical Estimates of $\\beta$ from the COVID-19 Trends and Impact Survey"),
subtitle = "Density Estimate of Distribution and Median by State" ) +
theme_c(axis.title=element_text(size=12),
axis.text.x=element_text(size=11)) ctis %>%
ggplot(aes(x = beta, fill=state)) +
geom_density(alpha=.6) +
scale_fill_viridis(option="mako", discrete=TRUE) +
theme_c()ctis %>%
group_by(state) %>%
summarize(mean=mean(beta,na.rm=TRUE)) %>%
pull(mean) %>%
mean()## [1] 0.2160426
ctis %>%
group_by(state) %>%
summarize(mean=mean(smoothed_wcli,na.rm=TRUE)) %>%
pull(mean) %>%
mean()## [1] 0.01602321
Compare Priors
source(here("R/05-state-analysis.R"))
beta_shape = get_shape(ctis_smoothed,
option="beta",
quiet=TRUE)
tibble(x = seq(0, 1, length = 10^(4)),
original = with(prior_params, beta_density(x,
mean = beta_mean,
sd = beta_sd,
bounds = beta_bounds)),
ctis= dbeta(x, shape1=beta_shape[1],
shape2=beta_shape[2])) %>%
pivot_longer(c(original,ctis))%>%
ggplot(aes(x=x, y = value, fill=name)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0, ymax=value,fill=name),
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $\\beta$"),
fill ='')s_untested_shape = get_shape(ctis_smoothed,
option="s_untested",
quiet=TRUE)
tibble(x = seq(0, 1, length = 10^(4)),
original = with(prior_params, beta_density(x,
mean = s_untested_mean,
sd = s_untested_sd,
bounds = s_untested_bounds)),
ctis= dbeta(x, shape1=s_untested_shape[1],
shape2=s_untested_shape[2])) %>%
pivot_longer(c(original,ctis))%>%
ggplot(aes(x=x, y = value, fill=name)) +
geom_line(alpha = .8) +
geom_ribbon(aes(x=x,ymin=0, ymax=value,fill=name),
alpha=.6) +
theme_c(legend.text=element_text(size = 10)) +
viridis::scale_fill_viridis(discrete=TRUE,
option = "rocket", begin=.3,end=.8) +
labs(x = "Value",
y = "Probability Density",
title = latex2exp::TeX("Definition of Prior for $Pr(S_1|untested)$"),
fill ='')